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ly-j ■ Abstract. A procedure for obtaining a "minimal" discretization of a partial 

differential equation, preserving all of its Lie point symmetries is presented. "Minimal" 
in this case means that the differential equation is replaced by a partial difference 
scheme involving N difference equations, where iV is the number of independent and 
\^ ' dependent variables. We restrict to one scalar function of two independent variables. 

As examples, invariant discretizations of the heat, Burgers and Korteweg-de Vries 
equations are presented. Some exact solutions of the discrete schemes are obtained. 
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1. Introduction 



The dynamics of most physical processes is described by differential equations. Lie 
group theory provides powerful tools for obtaining exact analytical solutions of such 
equations, specially the most fundamental ones that are usually nonlinear (the Einstein 
equations, the Yang-Mills equations, the Navier-Stokes equations, . . .). 

It is however quite possible that our world is actually discrete, at least at the 
microscopic level. If this is so then the fundamental equations are finite difference 
ones and the differential equations obtained in the continuous limit are actually 
approximations. Even if this is not so and no fundamental minimal space and time 
intervals exist, difference equations still play an important role in physics. On one 
hand, many physical phenomena are inherently discrete, such as vibrations in molecular 
and atomic chains, or phenomena in crystals. On the other hand, numerical methods 
for solving differential equations involve their discretization: the differential equations 
are replaced by difference ones, written on a lattice and these are then solved. 

In many cases, the symmetries of differential equations are more important and 
better known than the equations themselves since these symmetries reflect fundamental 
physical laws. Thus, when discretizing, or otherwise modifying dynamical equations, it 
is of interest to preserve the original symmetries. 

X Present address: Department of Mathematics, University of Minnesota, 526 Vincent Hall, 206 Church 
St. S.E., Minneapolis, MN 55455, USA. E-mail: valiq001@math.umn.edu. 
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This article is part of a general program that could be called "continuous symmetries 
of discrete equations" and that has been vigorously developed for the last 15 years, or so. 
Its overall aim is to turn Lie group theory into a tool for solving difference equations, just 
as it is for differential ones [1-24]. For recent reviews with extensive lists of references, 
see [18,24]. 

Different types of problems are treated in this program. A difference equation and a 
lattice may be given and the aim then is to solve the equations, using point symmetries, 
or generalized symmetries, as the case may be. Alternatively, as in this article, a 
differential equation may be given and the aim is to discretize it while preserving its 
important qualitative features, such as point symmetries. 

A formalism that unifies different approaches to symmetries of difference systems 
is one in which the difference equation and the lattice are described by a system of 
equations. These relate the dependent and independent variables evaluated in different 
points. This has been particularly fruitful for the discretization of ordinary differential 
equations (ODEs) [8,9,21]. 

The purpose of this article is to extend this approach to the case of partial 
differential equations (PDEs). We wish to approximate them by partial differential 
schemes (PASs) allowing the same Lie point symmetry group as the original PDEs. 
The schemes will be compared with related ones already existing in the literature [1- 
7,10,15]. 

A difference scheme of order K for an ODE consists of two equations relating K + 1 
points Xk and K + 1 values Uk (the discrete approximation of the function u(x) in Xk, 
evaluated in these points): 

E a ({x n+ j, u n+ j}j e j) = 0, a = 1,2, 

./ {M, M • 1 V - 1. .V}. n,M,NeZ, N > M. 

The points Xk are distributed along a line (the x-axis) and the two equations (1) must be 
such that if (x n+ M, • • • , x n +N-i, u n+ M, • • • , u n +N-i) are given, it is possible to calculate 

\ x n+N, u h+n}- 

The solution of the system will have the form 

x n = x(n, ci, . . . ,c 2 k), (2a) 
u n = u(n, ci, . . . ,c 2 k) (26) 

where ci, . . . , c 2 k are integration constants. 

In the continuous limit, one of equations (1) goes to an ODE, the other to an 
identity (like = 0). 

An essential observation is that in this approach the actual lattice (2 a) is not a 
priori given, but is obtained by solving the system (1). 

Point transformations for the system (1) will have exactly the same form as for the 
ODE, namely 

x — A g (x,u), u — Q g (x,u), (3) 
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where g represents the group parameters. They are generated by a Lie algebra of vector 
fields of the form 

V = Z(x,u)d x + <f>(x,u)d u . (4) 

The prolongation of the vector field V to all points figuring in the scheme (1) has the 
form 

n+N 

prV= ^ t( x k,u k )d Xk + 4>{x k ,u k )d Uk (5) 

k=n+M 

and the fact that the corresponding group transformations (3) will take solutions into 
solutions is assured by imposing 

prV[£ a ] =0, a = 1,2. (6) 

E 1= E 2 =0 

This approach has proven to be fruitful for large classes of ordinary difference 
schemes (OASs) [8,9,14,21]. It provides exact discretizations of first order ODEs [21] 
(i.e. OASs that have exactly the same solutions as their continuous limits) and second 
order OASs that can be exactly solved [8,9]. As pointed out earlier [15, 18,24] the use 
of point symmetries on fixed, nontransforming lattices is much less fruitful. 

In Section 2 we present the general symmetry preserving discretization of a scalar 
PDE with two independent variables. In spirit the method is the same as used for 
ODEs [8,9] and it leads to a system of 3 difference equations (rather than 2 as for 
ODEs). Sections 3, 4 and 5 are devoted to examples. The linear heat equation is 
treated in Section 3, the Burgers and Korteweg-de Vries equations in Sections 4 and 5, 
respectively. The final Section 6 is devoted to some conclusions and the future outlook. 



2. Invariant discretization of a partial differential equation 

For simplicity of notation we restrict ourselves to a PDE involving one scalar function 
of two variables u(x, t). It will be approximated by a difference equation on a symmetry 
adapted lattice. The lattice consists of points distributed in a plane. We will label these 
points by an ordered pair of integers P m ,n- We also introduce continuous coordinates in 
the plane and call them (x,t), though they do not necessarily correspond to space and 
time and are not necessarily cartesian coordinates. The coordinates of the point P m ,n 
will be 

(x m ,n,t m ,n), (m,n)eZ 2 , (7) 

see figure 1. 

The actual partial difference scheme will be a set of relations between the variables 
(x m ,n, t m ,n, u m ,n) evaluated at a finite number of points on the lattice. The first question 
is: how many relations and how many points do we need? 

We start from a given PDE 

E(x,t,u (n \x,t)) = 0, (8) 
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Figure 1. The (x,t) coordinates of points on a two-dimensional lattice. 



(10) 



where u^ n \x,t) denotes all the partial derivatives of u(x.t) up to order n. We assume 
that equation (8) is invariant under a group G of local Lie point transformations with 
a Lie algebra L realized by vector fields of the form 

V = t, u)d x + r](x, t, u)d t + <f>(x, t, u)d u . (9) 

We wish to approximate the PDE (8) by a system of finite difference equations 

Ea({Xrn+ji,n+j2itrn+ji,n+j2i U m+ji,n+j2}(ji,j2)£j) ~ 0> 

1 < a < N, where {0} C J C Z 2 , 

relating the quantities (x, t, u) at a finite number of points and invariant under the same 
group G as (8). The minimal number of equations needed in (10) is N = 3, determining 
the values of u, x and t at different points. If only 3 equations are imposed, then the 
solution of (10) will depend on a certain number of arbitrary functions of one variable 
(to, n or some combination of to and n). How many such functions will be present 
depends on the order of the system, i.e. the number of points corresponding to the set 
J. This in turn depends on the order of the PDE (8) and the precision we are looking 
for. 

It may be convenient to impose more than 3 equations, to specify the lattice 
to a larger degree. For example, in reference [15] the number of equations imposed 
was 5. Four of them specified the mesh and made it possible to move along the 
lattice in all directions from any point. The additional equations (N > 3) cause the 
system (10) to be overdetermined, so some compatibility conditions must be satisfied. 
These additional equations, if imposed, will play the same role as initial conditions, or 
boundary conditions for PDEs: they will partially, or completely specify the arbitrary 
functions involved. 

To illustrate the point, consider the linear PDE 

u xt = 0. (11) 
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We approximate it by the PAS [18] 

,-, 1 / ^m+l,n+l u m+l,n ^m,n+l u m,n \ n 

1 = 7 T7 — I Z Z j = u ' 

(12) 



-^2 ^m,n+l ^m,n 

This system can be solved explicitly to obtain 

t m , n = a(m), x mtn = P(n), (13a) 

Um,n = f(Xm,n) + g(t m ,n)- (136) 

Thus the general solution involves 4 arbitrary functions of 1 variable each. The lattice 
is orthogonal, the spacing of points in the two directions is arbitrary. We mention that 
the PAS (12) is invariant under an arbitrary reparametrization of t and x (i.e. of m 
and n), corresponding to the invariance of the PDE (11) under the infinite-dimensional 
group of conformal transformations. If we wish to impose that the lattice be uniform 
in both directions, we can add two further equations 

E& tm+l,n 2t mn + t m —\ :n 0, (14) 
E§ X m ^ n +\ 1x m ^ n -\- Xm^n—i 0. 

The solution (13 a) then restricts to 

t min = rm + to, x min = hn + x , (15) 

with Um^n as in (136) (to, t, xq, h are constants). The PAS given by (12) and (14) is no 
longer invariant under the infinite-dimensional conformal group. The symmetry group 
of the lattice is reduced to dilations and translations of t and x. 

Since we wish to obtain a symmetry preserving discretization of the PDE, the PAS 
(10) must be constructed out of invariants and invariant manifolds of the symmetry 
group G. The procedure is standard [1-10, 18,24] and we describe it only briefly. 

1. We specify a stencil i.e. choose the number and positions of the points to be used 
in the PAS. 

2. We prolong the vector field (9) to all points of the stencil (i.e. all the points figuring 
in the PAS) 

P r V = ^] {.£,m+j 1 ,n+j 2 dx m+ j 1>n+ j 2 + 1 lm+j 1 ,n+j 2 dt m+jl>n+J2 + 4 > m+j 1 ,n+j 2 9u rn+jlin+j2 } , (16) 
J 

where t, m+jun+h = £(x m+jun+j2 ,t m+jun+h ,u m+jun+h ) and similarly for i] and <p. 

3. We find the elementary invariants of G by solving the system of first order PDEs 

pr \[I({x m+ j lt n + j 2 , t m+ j ltTl+ j 2 , U m+ j ltn+ j. 2 \^j 1 ^ < zj)\ = 0, (17) 

where V is a general element of the symmetry algebra L of the PDE (8) (but the 
prolongation is the "discrete" one (16)). The algebra may be finite, or infinite 
dimensional. If we have dim L = I < oo then we choose a convenient basis V, 



O! 



a = 1,2, ... , / and (17) reduces to a system of / linear first order PDEs.) Using the 
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method of characteristics, we obtain a set of elementary invariants Ji, ...,J M . Their 
number is given by the formula 

H = dim M - rank Z, (18) 
where M is the manifold that G acts on, i.e. 

M ~ {{x m +j 1)n +j 2 , t m +j ljTl +j 2 , u m +j ljTl +j 2 y (j 1 ,j 2 )ej}- (19) 
Thus dim M = N x #J, where # J denotes the order of the set J and Z is the matrix 

^ {Cn+ji,n+j 2 > r lrn+j 1 ,n+j 2 i ( f } m+j 1 ,n+j 2 } (h ,h)eJ \ 

z = : ; (20) 

\ {Cn+ji,n+j 2 > r lm+j 1 ,n+j 2 i ( t ) m+j 1 ,n+j 2 ^ {h ,h)eJ / 

formed with the coefficients of the prolonged symmetry generators V a generating the 
basis of the finite dimensional Lie algebra. 

Since the quantities I±, . . . , 1^ form a basis of elementary invariants, any difference 
equation 

=0 (21) 

will be invariant under the group G. The equation (21) obtained in this manner is said 
to be strongly invariant and satisfies pr V a [_E] = 0, a = 1, . . . , /, identically. 

Other invariant equations can be obtained if the rank of the matrix 
Z is not maximal on some manifolds described by equations of the form 

E{\ x rn+j 1 ,n+j2>tm+j 1 ,n+j2i U m+j 1 ,n+j2}(ji,j2)£j) = that Satisfy 

prV a [E] =0, a = l,...,/. (22) 

Such equations are said to be weakly invariant. In practice we can start by computing 
the invariant manifolds and this can facilitate the computation of the set of fundamental 
invariants. 

In this article we will always choose the minimal number of equations needed, i.e. 
N = 3 in (10). By construction these equations all satisfy 

prV[E a ] =0, a =1,2, 3, (23) 

E 1 =E 2 =E :i =0 

i.e. they are invariant under the group G. 

When constructing the PAS out of the invariants Ii,...,I^ and the invariant 
manifolds some choices must be made. The most important constraint is given by 
the continuous limit. We impose that E\ — reduces to the PDE (8) and E 2 — E 3 — 
reduce to identities. Further constraints may come from the boundaries, or initial 
conditions imposed on solutions of the PDE that we are solving, or from the precision 
that we desire. 

The calculation of group invariants described above is completely analogous to the 
one used in the continuous case [25,26]. The Lie algebra approach that we use can be 
replaced by the Lie group one using moving frames [27, 28]. 
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If the symmetry group of the original PDE is infinite-dimensional, some 
modifications of the procedure are required. In particular, if the PDE is linear, then 
an infinite-dimensional pseudogroup corresponding to the linear superposition principle 
is always present. In this case we can restrict ourselves to invariants of the finite 
dimensional subgroup of the symmetry group and then require that the PAS formed 
out of the invariants be linear in u. 

3. The linear heat equations 

The linear heat equation 

iH = u xx (24) 

is a much used example when Lie group theory is applied to the study of differential 
equations. Its symmetry group was already known to Sophus Lie and is reproduced in 
virtually every book on the subject (see e.g. [25]). A basis for its symmetry algebra 
consists of 

Vi = d x , V 2 = dt, V 3 = ud u , V 4 = xd x + 2td t , 

(25a) 

V 5 = 2td x - xud u , V 6 = Atxd x + At 2 d t - (x 2 + 2t)ud u , 

\ a = a(x,t)d u where a t = a xx . (256) 

Like everybody else, we use the heat equation (24) as an example, precisely because 
it is linear, as reflected in the infinite dimensional algebra (256) and has a large and 
interesting finite dimensional subalgebra (25a) of the symmetry algebra L. 

Our aim is to discretize (24) while preserving the entire Lie point symmetry algebra 
(25). Several similar discretizations exist in the literature. In [1] and [7] the authors 
first introduce Lagrangian variables and represent the equation (24) as a system of two 
equations. These are then discretized preserving all the symmetries (25 a), not however 
linearity. A different point of view was taken in [11,13,16]. Equation (24) (and any 
other linear equation) is discretized on a uniform orthogonal nontransforming lattice. 
Linearity is preserved, but the Lie algebra of point symmetries (25) is replaced by an 
isomorphic Lie algebra of difference operators. Thus, the Lie point character of the 
transformations is given up while all the algebraic consequences of Lie symmetry are 
preserved. Here we shall directly apply the procedure outlined in Section 2. 

3.1. Invariant schemes for the linear heat equation 

Before computing a set of elementary discrete invariants, let us introduce a convenient 
notation for points on the lattice, to be used for the heat equation and all further 
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(x—,t) (x-,t) (x,t) (x+,t) (x++,t) 

Figure 2. Lattice for the invariant discretization of the heat and Burgers equations 
(6 points) and the KdV equation (10 points). 



examples in this article. We put: 

{%m,ri) tm^n, Wm,n) = (x, t,"u), 
{%m+l,n-> tm+l,ni ^m+l,n) = {%it,U), 

(%m, rail i ^m, nil i ^m, nil) = (%±i t±i U±) , / rir ,\ 

(zo) 

(^m,n±2) ^m,ni2, U m ,n±2) = (^ii, *±± , U±±) , 
(^m+l,nil) ^m+l,n±l) u m+l,n±l) = (&±it±,U±), 
(Xm+l,n±2, *m+l,n±2, ^m+l,rai2) = (£±±, i±±, W±±), 

and introduce the steps 

/i± = ±(x± — x), fo±± = ±(x±± — x±), /i± = ±(x± — x), 

fo±± = ±(x±± — x±), er = x — x, ct + = x + — x + , (27) 

T = i-t, T± = ±(t±-t). 

(this notation is standard in numerical analysis). 
It is easy to see that the equation 

T + = (28) 

is invariant under the entire group G generated by the algebra L of (25) (i.e. we have 
pr V[T + ]| T+=0 = for all V e L). We will include (28) in all of our invariant schemes 
which means that we will always have horizontal time layers. Equation (28) implies that 
in (26) we have 

t_| |_ = = t = t— = t , t_| [_ = t^ = t = t_ = t ; 

as indicated on figure 2. The fact that the time layers are horizontal is particularly 
convenient in numerical simulations. 

For the heat equation we will not need all 10 points of figure 2. We will use 6 of 
them and in view of (28) we can restrict to the 14 dimensional space with coordinates 

(x_, X, X + , X_, X, X + , t, t, M_, U, U + , U_, u, u + ). (29) 
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Imposing invariance under the group generated by the 6 dimensional algebra (25 a) we 
obtain 8 elementary invariants 



K 










h 


h- 


1 




u+ 




rh + 


u 


exp 


At 


U+ 




\K 




exp 


At 


u 





h+h A 



h = — exp 

u 

1 8 = — exp 

u 



a 



At 



<TU 

exp 

h + u 

^(2a + h_) 



(30) 



This set is equivalent to the one used in references [1,7], but we find the set (30) more 
convenient when imposing linearity. 

Using the invariants (30) we obtain an explicit invariant scheme that is linear in u 
by setting 

A 



7-3/2 j j 



' 3 J 4 - -13 - (-^5 + h) exp 
In terms of the variables (x, t, u) we have 
1 



-2, 



0, 



h 



(31) 







ff>. 


' a ' 


/ —u exp 
/ h 


At\ 




1 




hh 




T+ 




h + 



.)- 



^-{2a -h + h) 
At 



2u + exp 



^-(2a + h-h) 
At 



u + exp 
= 0, 

= h- = h. 

An invariant implicit scheme, also linear in u, is obtained by setting 

V 



,1/2.-1 
'3 1 A 



(I 7 + 1 8 ) exp 



2, 



T + = 0, 



1, 



(32 a) 

(326) 
(32 c) 

(33) 



which in terms of the discrete variables gives 



1 



u — 



-u exp 



a 

Arl 



u + exp 



h 



At 



[2a - h + h) 



— 2u + u^ exp 



h 



(2a + h-h) 

At 



(34) 



hh^ 
T + = 0, 

h + = h- = h. 

We recall that in this context "explicit" means that one value of u on a higher 
time level is express in terms of values of -a at a previous time t. "Implicit" means that 
several values of u at time t are calculated simultaneously (for x+, x and £_). 

Notice that unlike for the standard explicit and implicit discretization 
u — u u + — 2u + U- u — u u + — 2u + -u_ 

— = — s> — ■ — = — p — ' (35) 

it is not possible in the invariant case to have schemes involving just h or h. Indeed, in 
the invariant schemes both h and h are present. This is a consequence of the use of J 3 
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when generating the schemes. In the particular case when a is chosen to be zero, we get 
h = h and then the invariant schemes reduce to the standard discretizations (35) on an 
orthogonal lattice. It is important to realise that the choice a = is not an invariant 
one. In fact, it is not invariant under the transformations generated by V5 and Yq. 

The two equations for the lattice of the invariant schemes (32) and (34) are easily 
solved and give 

t m , n = 7H, ^m,n = h(m)n + x (m), (36) 
where 7(777), h{m) and x (m) are arbitrary functions. 

3.2. Continous limit of the invariant schemes 

We now compute the continuous limit of the explicit invariant scheme (32) to first of all 
obtain an important condition on the limit of the ratio a/r and show that the scheme 
then coverges to the heat equation. We will not present the calculations for the implicit 
schemes since they are quite similar. 

Clearly, the equations for the lattice (32b) and (32c) go to 0=0 in the continuous 
limit. We must show that (32a) goes to (24). To compute this limit, we introduce 
infinitesimal parameters that go to zero in the continous limit. Since the steps induced 
by the discrete variable n must go to zero independently of those induced by the discrete 
variable 777, we must have 

7(777) = €7(771), h{m) = 5h(m), (37) 

where e and 5 are independent infinitesimal parameters and 7(777) and h(m) are finite 
quantities. By the definition of x m , n in (36) we have 

a = (h(m + 1) — h(m))n + x (777 + 1) — 2:0(777) 
= 5(h(m + 1) — h{m))n + x (m + 1) — £0(777). 
Since, a must go to zero when r goes to zero, a must converge to zero as 

a = 8e k ai (777)77 + e l a 2 (m), (38) 
with k, I > 0. With these considerations we have 

in the continous limit. By writing 

h = h + a + -a, (39) 
the expansion of 




h ' 

in a Taylor series is well defined. By developing the left hand side of (32 a) we find 

+ -r^u+-u x + u t + 0(-(^—) ,—,r,a). (40) 



2 rh 4r 2 r r V h J r 
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Equation (40) written in terms of the parameters e and 5 gives 

l -e k ~ l + J(5e fc -Vi(m)n + e^V 2 (m)) 2 M 

+ (5e fc -V!(m)n + e^V 2 (m)) Ma; + M< + 0(e 2fc - 1 , e 2 '- 1 ) 
For the coefficient in front of in (41) not to diverge in the limit we must impose 

k, I > 1 (42) 
For the invariant scheme to converge we must have 

lim - < oo. (43) 

The expansion of the right hand side of (32 a) in a Taylor series gives 

1 <7_i_ — a a 2 a . 9l . . . 

+ —2 u + -u x + u xx + 0(e 2 k-l,5). (44) 

Equating (40) to (44) we get 

u t = u xx + 0(e 2k -\e 2l -\6), 
which converges to the heat equation when e and 5 go to zero. 

3. 3. Exact solutions of the invariant schemes 

One of the interesting aspects of generating symmetry-preserving schemes is that group 
theory can be used to obtain nontrivial exact solutions, using the fundamental concept 
that the symmetry group maps solutions to solutions. 

One obvious solution of the schemes (32) and (34) is the linear solution 

u m ,n = ax min + b, (45 a) 

defined on the orthogonal lattice 

t m ,n = l(m), x m , n = hn + x , (45b) 

where h and x are constants. Indeed, as mentioned earlier, on such a lattice a = so 
the invariant discretizations go to the standard discretizations (35) and it is well known 
that (45a) is an exact solution. If we want, we can fix / ~f(m) to be equal to rm + to, 
where r and t are constants, so that the lattice is rectangular. 

Given a discrete solution (x,t,u(x,t)) new solutions can be obtained by acting on 
the known solution with the symmetry group generated by (25 a). Hence, 

e ei (x + e 1 ) +2e 5 e 2e4 (t + e 2 ) ~ e 2e4 (t + e 2 ) 



x 



l-4e 6 e 2 ^(t + e 2 ) ' 1 - 4e 6 e 2 ^ + e 2 ) ' 



~(~7\ 1 ( - cs x-2e 6 t e-^t 
u(x,t) — — u e b ~, - — e 2 exp 

VT+4eJ I 1 + 4eet 1 + 4eet 



e 5 x — e\t + e 6 x z 

e 3 ~ 

1 + 4e 6 t 



(46) 
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is also a solution of the invariant schemes (32) and (34). If the general transformation 
is applied to the linear solution, we conclude that 

e £4 (/m + x + ei) + 2e 5 e 2e4 7(m) 
X l-4e 6 eM 7 (m)+e 2 ) ' 1 j 

T e 2e4 ( 7 (m) + e 2 ) ? 



4e 6 e 2e4 (7(m) + e 2 ) ' 

e 5 x — e\t + e 6 :r 2 



~l~ 7\ 1 C -65 X ~ ^ , M 

t) = — _ (ae 5 ~ + oj exp 

Vl + 4e 6 t l+4e 6 t 



£3 



1 + 4e 6 t 



,(47c) 



is a solution. 

Furthermore, the limit condition (43) for the ratio a/r is preserved under the 
general symmetry transformation (47a) and (476), since 

a = x-x^ = a^ e4 _ ^ + + ^ + 4ege£4(x + £i) (4g) 
T t-t T 

which does not diverge in the continuous limit by hypothesis on the ratio a/r. 

The fundamental solution of the heat equation can be obtained by successively 
applying the symmetry transformations exp[eV 6 ], exp [ln( v /0rt 2 )V 3 ] and exp[l/(4e)V 2 ] 
to the constant solution u = b ^ and on the lattice x = hn + xo, t = 7(771) we get 

(49 a) 



// = \ I — - exp 

4nt 



—x 



At 



defined on the lattice 



x ~ I 

% = 7~, t= ^r-i 7-^- (496) 

l-4et' 4e(l-4et) v ; 

By taking e = 1/4 and 7(771) = \ (l — m ^ +tQ j , we get the lattice 

x = (hn + x )(rm + t ), t — rm + to, (50) 

where h, xo, r and to are constants. The evolution of the lattice is shown in figure 3. 

Another interesting solution that can be obtained is the traveling wave solution. 
This one is obtained by applying the symmetry transformation exp[cV 5 ] 

u = (a(x — let) + 6) exp[— cx + c 2 t], (51a) 

and is defined on the lattice 

x = hn + xo + 2ct, t — 7(771). (516) 

By taking 7(771) = rm + t , we obtain a lattice in which the points in x move linearly 
as a function of time. 
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Figure 3. Lattice for the fundamental solution of the heat equation, h = 0.05, xq = 0, 
t = 0.005, t Q = 10. 



4. Burgers' equation 

In this section we shall derive invariant schemes for the Burgers' equation 

v t + vv x = v xx . 



(52) 



This equation is well known to be related to the heat equation (24) by the Cole-Hopf 
transformation 



v(x, t) 



u(x, t) 



(53) 



Since the transformation (53) is not a point one, the symmetry algebras of the Burgers' 
equation and the heat equation are not isomorphic. Indeed, the symmetry algebra of 
(52) is five dimensional and is spanned by the vector fields [25] 

Vi = d x , V 2 = d t , V 3 = td x + d v , 



V 4 = xd x + 2td t — vd V) V 5 = txd x + t 2 d t + (x — tv)d v . 



(54) 



4-1. Invariant schemes 



Before computing a set of fundamental discrete invariants let us mention the fact that 
the equation 

T + = (55) 



Discretization of partial differential equations preserving their physical symmetries 14 



is again weakly invariant under the group corresponding to (54). The set of elementary 
invariants involving the same discrete points as for the heat equation (29) is 



h + h + T h + h^ 



h- h- t 

h = h + h-(v+ -v~), h = h+h-{v+ -v~), 

h = h +(^- v )i h = h+ - vj , 



(56) 



where 

, v + — V _ v — V _ A , v + — V ^_ v — v_ 

v x = —, ) V x = —, > V x = — > v x = — • ( 57 ) 

h+ h- h + h- 

From the set of invariants (56) and the weakly invariant equation (55) we derive 
invariant explicit and implicit schemes. The explicit scheme is obtained by putting 

{2h-I 7 )h-hh = h, T+ = 0, h = l, (58) 

or in terms of the discrete variables 

T+ = 0, (596) 

h + = h_ = h. (59c) 
The implicit scheme is obtained by putting 

(J 6 - 2I 7 )I 3 - hh = h, T + = 0, I 2 = l, (60) 



T + = 0, (616) 

h + = h^ = h. (61c) 

As for the heat equation, the schemes generated involve h and h. Furthermore, the 
computation of the continuous limit of the invariant schemes (59) and (61) gives the 
same condition on the ratio d/r as for the heat equation, namely (43). To see this, we 
compute the continuous limit of the explicit scheme (59). The calculation is similar for 
the implicit scheme (61). As for the heat equation, we suppose that the steps in t and 
x are given by (37) and (38), where we immediately assume that k, I > 1 in (38) . First 
of all, it is clear that in the limit the equations specifying the lattice (596) and (59c) 
go to the identity = 0. The development of (59a) in a Taylor series in terms of the 
infinitesimal parameters e and 8 gives 



(7 a I a + — a \ v a I a + — a\ I a + — a 

T 2 T\hT/TT Xt \ hr J \ hr 

o~ ( o+ — a\ a v a . 2k _, 27 _i . 

vv x --v x + 7 = v xx + 0(e 2k \e 21 \S). 



(62) 



t \ hr j T T T 
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In (62) we have explicitely written all the terms that do not converge to zero in the 
limit, but after canceling the superfluous terms and making e and 5 go to zero we 
retrieve Burgers' equation (52). 

In the particular case where the steps in t are constant and a is zero, the invariant 
discretization reduces to standard one, namely 

v — v , v+ — v~ 

r x h 

t = rm + t , x = hx + x , 

for the explicit scheme and 

v — v _ , vt 



+ vvZ = 1 



r x h 1 

t = rm + t , x = hx + x , 

for the implicit one. 

The invariant schemes just generated are closely related to those obtained by 
Dorodnitsyn and Kozlov [6]. The invariant schemes (59) and (61) are not uniquely 
defined, the discrete approximations of the Burgers' equation (59 a) and (61a) can be 
defined on other lattices. For the explicit scheme we can replace (59 c) by the invariant 
equation 

h = 0. (63) 
On this new lattice, (59a) becomes 

-hh = h- (64) 
In terms of the variables (x, t, v) we get the invariant explicit scheme 

v — v h- . , , 
t hi 

(65) 

a = tv. 

For the implicit symmetry-preserving scheme we can replace (61c) by Ij = to obtain a 
scheme similar to those in [6] . In the case of the Burgers' equation the invariant schemes 
in [6] can be seen as particular cases of the schemes we obtained. The major difference 
between the invariant schemes in [6] and (59) and (61) is that the steps in x do not to 
have to be uniform in general. 

4-2. Exact solutions 

Similarly as for the heat equation, it is clear that the constant solution 

v(x, t) = vq, (66a) 
defined on the orthogonal lattice 

t m ,n = l(m), x m<n = hn + x , (666) 
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is an exact solution of the invariant schemes (59) and (61) (t> , h and x are constants). 
By choosing 7(771) = Tin + t we get a standard rectangular lattice. Since the 
symmetry group of the Burgers' equation is not very rich, just one new solution can be 
obtained from the constant one. By applying successively the symmetry transformations 
exp[eVs], exp[t>o/eVi] and exp[l/eV 2 ] to (66a) and (666) we obtain 

v = x/t (67a) 

definied on the lattice 

~ _ ex m , n + (1 - et m , n )v ~ _ 1 

X m ,n (1 + \ ' L m,n (1 + \ ^ ' 

If we choose vq = 0, take e = 1 and pose 7(771) = 1 — l/(rm + 1) in (666), we conclude 
that the solution (67a) is an exact solution of the invariant schemes on the lattice 

t m ,n = Tm + to, X min = (hn + Xo){T17l + t ). (68) 

The same lattice over which the fundamental solution of the heat equation is exact. 
4-3. Burgers' equation in potential form 

The invariant discretization of the Burgers' equation in potential form 

w t H y~ = Wxx ( 69 ) 

is readily obtained from that of the linear heat equation. Indeed, since the PDE (69) is 
related to the heat equation (24) by the point transformation 

u = exp |-— J , (70) 

the results of Section 3 can directly be used to generate invariant schemes and exact 
solutions via the transformation (70). For instance, the symmetry algebra of (69) is 
spanned by 

Vi = d x , V 2 = d t , V 3 = td x + xd w , V 4 = xd x + 2td t , 

x 2 ( 71a ) 
V 5 = txd x + t 2 d t + {— + t)d w , V 6 = d w , 

[w~\ 

V a = a(x,t)exp y— J d w , (716) 
where a is a solution of the heat equation, a t = a xx . 



4-3.1. Invariant schemes Since the transformation (70) does not affect the independant 
variables, it is clear that the equation (28) is still a weakly invariant equation. We can 
directly use the set of elementary invariants (30) and the transformation (70) to obtain 
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a set of elementary invariants on a flat time layer scheme [23] 



h 

h 



h+ 

exp 
exp 



h = 



h+ 
h_ 



h = 



h+h^ 



h 



exp 



+ 



w + 



— w h + , n , . 



w + —w h+ . f . 



exp 



exp 



2 

u> — u). 



w — w a 
2~ + 47 

- — (2a + h. 
At 



(72) 



By using the same invariant expressions (31) and (33) we obtain invariant schemes for 
(69). An explicit scheme is 




w 


a 2 ' 




" w 




~1> 


+ 47 


— exp 


~2. 





-^ + -(2a-h + h) 



—2 exp 



w 
2 



+ exp 



W- h , n , fs 
— --(2a + h-h) 



(73) 



T + = 0, 
h + = h- = /i. 
An implicit one is given by 



2" 




a 
4r~ 



^— exp 
hh 



w + h 



--± + -(2o-h + h) 



2 exp + exp 



~-^ + h-h) 



(74) 



T + = 0, 

h + = h_ = h. 

The invariant schemes (73) and (74) converge to the differential equation (69), if 
the condition (43) on the lattice is satisfied. In contrast to the invariant schemes for 
the heat equation (32) and (34), which give the standard discretization when a = 0; 
the invariant schemes (73) and (74) give unusual schemes on a rectangular lattice. For 
example, when a = the explicit scheme (73) becomes on a rectangular lattice 

exp [-1 ] - exp [-f ] exp [-^} - 2 exp [-f ] + exp [-^} 

r h 2 ' (75) 

t = rm + to, x = hn + Xq, 

which is quit different from the standard discretization 

2 



W 



— W 1 

— +2 



w 



w + — 2w + w. 



h J h 2 

Finally, let us mention that the exact solutions obtained in Section 3.3 are 
mapped by (70) to exact solutions of the invariant schemes (73) and (74). Since the 
transformation does not affect the independant variables the solutions obtained via (70) 
are exact on the same lattice. 
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5. Korteweg-de Vries equation 

As a last example, we consider the Korteweg-de Vries equation 

u t = uu x + u xxx . (76) 
This equation has a four dimensional symmetry algebra spanned by 
Vi = d x , V 2 = dt, Vs = td x -d u: V 4 = xd x + Md t -2ud u . (77) 

5.1. Invariant schemes 

Again, the equation (28) is weakly invariant and we can discretize (76) on a lattice 
with flat time layers. This time we make use of all points shown in figure 2 in order to 
approximate correctly the third order derivative in x. The usual computation gives the 
invariants 



I - h+ 


T K+ 
K 


h = 


h- 
h—' 


/ - h+ 


T K+ 

h+ 


h = 


h- 
h.J 


i - h+ 

h+ 


h 3 

T - + 
r 


h = 


a + tu 
h + 


ho = tu~~ ) 


hi = tu~, 


h2 = 


ru^ 


h 3 = TU+ + , 


Iu = h 2 + (u-u), 


h 5 = 


TU", 


he = ru~, 


I 17 = ru+, 


hs = 


TU++, 



(78) 



with 

Mil — Mi U_ — U U++ — U+ „ U_ — U .„ . 

«r = ^r — ^ u *~ = — - h — , K + = — ± > fl r = — £ — • (79) 

n ++ n — h ++ h 

An invariant explicit scheme approximating (76) is obtained by putting 

Iu = hh h2 ~ — + ^(/i3 - In - hi + ho), T + = 0, h = l. (80) 

In terms of the original variables we have 

u — u u + — U- u ++ — 2u + + 2w_ — u au + — 

— = U —2h— + 2tf + T—2h—> (81a) 

T + = 0, (816) 

h + = h^ = h. (81c) 

Furthermore, an implicit scheme is obtained by considering the combination of 
invariants 

J 14 = (J 9 + I^I 14 )I 8 Il7 ~ Iw + i(/ 18 - I 17 - h 6 + h 5 )I 7 , T+ = 0, h = 1, (82) 
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i.e. 

u — u — u_ u ++ — 2u< + 2m_ — u au + — u_ . „ . 

= u — — — - h — — 1 — , (83a) 

t 2h 2h 3 r 2h 

T + = 0, (836) 
h + = h^ = h. (83c) 

The invariant schemes include the standard discretizations. Indeed, if a equals zero 
then the last term on the right hand side of (81a) and (83a) vanishes, which is just the 
usual discretization of the Korteweg-de Vries equation. In order for the equations (81a) 
and (83a) to converge to (76) the condition 
,. o- 

hm — < oo 

t^O T 

must be imposed. 

As for the Burgers' equation, by considering the invariant schemes (81) and (83) 
on a lattice depending on the solution we recover similar schemes as those obtained in 
[4]. For the explicit scheme (81), if we replace (81c) by 

h = 0, (84) 
the equation (81a) then becomes 

/l4 = ^(/l3-/l2-/ll + /l0)- (85) 

In term of the variables (x, t, u) we have 

T h 2 + 

T + = 0, ( 86 ) 
a = —tu. 

For the implicit scheme (83) we just have to replace (83c) by I 9 + Iulg 1 = to get the 
invariant implicit scheme 

u-u = (ut + -ut) - (u--u--) 
r h 2 + 

(87) 

a = —tu. 



5.2. Exact solutions 

As with all previously considered invariant schemes, the constant solution defined on an 
orthogonal mesh is an exact solution. Another exact solution of the invariant schemes 
corresponds to the solution invariant under the Galilei transformation generated by the 
vector field V3. This exact solution is obtained in a non-obvious fashion. The idea is 
the following one. Given the partial differential equation that we approximate, we can 
perform a symmetry reduction to obtain an exact solution u(x,t). Given this solution, 
we can insert it in the discrete equation approximating the original differential equation, 
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which will give an equation relating the independant variables. If the equation obtained 
can be solved using the existing liberty in the lattice, then we have generated an exact 
solution of the PAS. 

In the continuous case solution invariant under the infinitesimal symmetry generator 
V 3 is 

u = - X -. (88) 
By replacing this solution in (81a), we get 

x = —. (89) 

T 

If we take the evolution in time to be 

t m ,n = rm + t , (90) 

the equation (89) can be solved for x and gives 

x m ,n = (hn + x )(rm + t ). (91) 

Hence, we conclude that (88) is an exact solution of the explicit scheme (81) on the 
lattice defined by the equations (90) and (91). 



6. CONCLUSIONS 



The discretization procedure presented in Section 2 and applied to specific examples in 
Section 3, 4 and 5 is a "minimal" one. By that we mean that the PDE (8) is replaced 
by a PAS (10) involving only N = 3 difference equations. The PAS is by construction 
invariant under the symmetry group of the original PDE. The general solution of the 
PAS (10) will depend on several arbitrary functions of one (discrete) variable. A specific 
solution u(x, t) and a specific lattice on which the solution is valid is obtained once these 
arbitrary functions are specified. 

The freedom inherent in these arbitrary functions can be used in different ways. For 
instance, we can require that certain physically important solutions of the PDE should 
also be exact solutions of the PAS. This was done in Section 3 with the fundamental 
solution of the heat equation. 

An alternative possibility is to use the invariants of the symmetry group G to add 
further invariant equations to the "minimal" PAS (10). This creates an overdetermined 
system of difference equations, restricting the freedom in the solutions. Since the system 
of N equations with iV > 4 is overdetermined, their compatibility must be assured. The 
idea of imposing an overdetermined system of invariant difference equations has already 
been explored earlier [3, 7, 15] with the aim of further restricting the form of the lattice. 

All PDEs studied in this article are evolution equations of the form 

Ut = f (x,t, u, u x , u xx , u xxx ) (92) 

and we made use of the fact that their symmetry groups always allowed horizontal time 
layers. 



Discretization of partial differential equations preserving their physical symmetries 21 

Other types of equations, requiring different lattices, will be studied in the future. 
The extension to more than one dependent variable and more than two independant 
ones is obvious. 

This research program, "continuous symmetries of discrete equations" has both 
physical and mathematical aspects that must be further pursued. One conclusion is that 
different discrete physical systems require different approaches. It was shown elsewhere 
[12, 13] that symmetries of linear theories can be studied in terms of commuting 
difference operators on fixed lattices. The results of this and related papers [1-10, 
14, 15, 17,21] indicate that at least for nonlinear discrete phenomena, the lattice should 
be considered as a dynamical one evolving together with the solution and described by 
an invariant system of difference equations. 

From the mathematical point of view we see as the greatest challenge the 
development of symmetry adapted numerical schemes that give better results than 
standard numerical methods. 
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